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Abstract 

We consider the problem of scheduling arrivals to a congestion system with 
a finite number of users having identical deterministic demand sizes. The con¬ 
gestion is of the processor sharing type in the sense that all users in the system 
at any given time are served simultaneously. However, in contrast to classical 
processor sharing congestion models, the processing slowdown is proportional to 
the number of users in the system at any time. That is, the rate of service ex¬ 
perienced by all users is linearly decreasing with the number of users. For each 
user there is an ideal departure time (due date). A centralized scheduling goal 
is then to select arrival times so as to minimize the total penalty due to devia¬ 
tions from ideal times weighted with sojourn times. Each deviation is assumed 
quadratic, or more generally convex. But due to the dynamics of the system, the 
scheduling objective function is non-convex. Specifically, the system objective 
function is a non-smooth piecewise convex function. Nevertheless, we are able to 
leverage the structure of the problem to derive an algorithm that finds the global 
optimum in a (large but) finite number of steps, each involving the solution of 
a constrained convex program. Further, we put forward several heuristics. The 
first is the traversal of neighbouring constrained convex programming problems, 
that is guaranteed to reach a local minimum of the centralized problem. This is a 
form of a “local search”, where we use the problem structure in a novel manner. 

The second is a one-coordinate “global search”, used in coordinate pivot itera¬ 
tion. We then merge these two heuristics into a unified “local-global” heuristic, 
and numerically illustrate the effectiveness of this heuristic. 
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1 Introduction 


Users of shared resources are frequently faced with the decision of when to use the 
resource with a view of trying to avoid rush hour effects. Broad examples include, 
workers taking their lunch break and attending a cafeteria; people entering and va¬ 
cating sporting events; and commuters using transportation networks. In many such 
situations the so called rush-hour game is played by all users acting individually. On 
the one hand, each user typically has an ideal arrival/departure time, while on the 
other hand, users often wish to avoid rush hour so as to minimise congestion costs. 
These general types of scenarios have received much attention through the transporta¬ 
tion community, [I], the queueing community (see [8] or p84 of ra for a review) and 
more specifically within the setting we consider in this paper [20]. 

While understanding social strategic (game) behaviour is important, a complemen¬ 
tary analysis is with regards to the social optimum (centralised scheduling decisions). 
These types of situations occur often in manufacturing, appointment scheduling, edu¬ 
cation and service. Most of the research on scheduling methodology does not consider 
processor sharing but rather focuses on the situation where resources are dedicated, 
see [18]. In this paper, we put forward a novel scheduling model, that offers a simple 
abstraction of a common scenario: Jobs may be scheduled simultaneously, yet slow 
each other down when sharing the resource. In this respect our model is related to 
the study of scheduling problems with batch processing, see [19]. However, from a 
mathematical perspective, our model, results and methods do not involve the classical 
discrete approaches but rather rely on piecewise affine dynamics with breakpoints. This 
type of behaviour resembles Separated Continuous Linear Programs, as in [53], and is 
often used to solve optimization problems associated with fluid multi-class queueing 
networks (cf. [2], [T6]). 

A standard way of modelling resource sharing phenomena, is the so-called processor 
sharing queue, see for example [9]. In such a model, given that at time t there are q(t) 
users in the system, the total fixed service capacity, f3 > 0, is allocated, such that each 
user receives an instantaneous service rate , 


d« ( ‘ ) ) = ity 


( 1 . 1 ) 


Such a model then captures the relationship of the arrival time of a user, a, the depar¬ 
ture time of a user, d and the service demand, £ through 





The aggregate throughput with q users in the system is the product qv(q). For the 
processor sharing model (II.ID . this is obviously the constant (3. However, in practice, 
the aggregate throughput is not necessarily constant with respect to q(t). In many 
situations, most notably in traffic and transportation scenarios, users inter-play in a 
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Figure 1: The service rate and aggregate throughput as a function of the number of 
users in the system. Parameter values: ft = 100 and a = 5 


complicated manner. In particular, in the classic Greenshield fluid model, (see for 
example HU or M) the aggregate throughput is not monotone in the number of 
users and even exhibits a traffic jam effect. The simplest model, describing such a 
phenomenon is 

v{q(t)) =fi-a(q(t) - l), (1.2) 

which is a discrete variation of Greenshield’s modefl. With a single user in the system, 
(II.21) yields the free flow rate ft which coincides with (II.IF Then for each additional 
user, there is a linear slowdown of a > 0 units in the rate. See Figure |T] for a simple 
illustration. Note that in road networks, much research has focused on the so-called 
fundamental diagram for networks, such as in [ 6 j. Indeed Figure |T]-b resembles a 
fundamental diagram. 

Our scheduling problem is to centrally choose arrival times a = (op,..., a n)' in an 
effective manner, where N is the number of users. In this paper we assume that all 
users share the same service demand, l. In our objective, user i incurs a cost of 

(di - d *) 2 +7 (di- Oi ), 

where di is his departure time and d* is the ideal departure time (due date) and 7 
captures tradeoff between meeting the due date and sojourn time costs. The total 
costs incurred by all users is then the sum of individual user costs. 

If there was no congestion (say due to d* being well separated), an ideal choice is 
dj = d* — (./ft. But in general, users interact, so the scheduling decision needs to take 
this interaction into account. If, for example, 7 = 0 and d* = d* for all i, then the 
problem is trivially solved with zero cost by setting 

,, t 

a i = d — 77 - TTj 

ft — ce(N — 

^^Note that in queueing theory, situations where v(-) is not as in (ED but is rather some other 
function are sometimes referred to as generalized processor sharing. See for example [5]. Generalized 
processor sharing has also taken other meanings over the years, so sometimes there is confusion about 
the term. 
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Here since sojourn time does not play a role, sending all users simultaneously will 
imply they arrive simultaneously after being served together at the slowest possible 
rate. Continuing with the case of 7 = 0, if now users do not have the same d*, then 
attaining zero costs is still possible. In fact, we show in the sequel, that in this specific 
case (7 = 0 ) the optimal schedule can be computed efficiently (in polynomial time). 

At the other extreme consider the case where minimising sojourn times is prioritised 
over minimisation of due dates (e.g. if fuel costs are extremely high). This corresponds 
to 7 ~ 00 . While for any finite 7 , it is possible that an optimal schedule allows overlap 
of users, an approximation for the case of large 7 is obtained by enforcing a schedule 
with no overlap ( q(t ) < 1 Vf). This is because overlaps have a very large sojourn time 
cost relative to the possible reduction in quadratic deviation from desired departure 
times. Now with such a constraint, the problem resembles a single machine scheduling 
problem with due date penalties. This problem has been heavily studied (see for 
example [3] or [ 2 Tj ). In our case, in which users have identical demand, finding the 
optimal schedule is a convex quadratic program and can thus be solved in polynomial 
time. We spell out the details in the sequel. 

Setting aside the extreme cases of 7 = 0 or 7 « 00 , the problem is more complicated. 
While we do not have an NP-hardness proof, we conjecture that finding the optimal a is 
a computationally challenging problem. In the current paper we handle this problem 
in several ways. First we show that departure times depend on arrival times in a 
piecewise affine manner. We find an efficient algorithm for calculating d*(a). We then 
show that the total cost is a piecewise convex quadratic function but generally not 
convex, i.e. there is a large (but finite) number of polytopes in where within each 
polytope, it is a convex quadratic function of a. This is a similar formulation to that of 
the piecewise-linear programming problem presented in [23j, which is known to be NP- 
liard. The structure of the total cost yields an exhaustive search scheduling algorithm 
which terminates in finite time. 

We then put forward heuristics. The first heuristic, which we refer to as the local 
search, operates by solving a sequence of neighbouring quadratic problems until find¬ 
ing a local minimum with respect to the global optimization. The second heuristic 
performs a global search over one coordinate (arrival time of a single user), keeping 
other coordinates fixed. This is done in a provably efficient manner. In particular, we 
bound the number of steps in each coordinate search by a polynomial. It then repeats 
over other coordinates, cycling over all coordinates until no effective improvement in 
the objective function is possible. In case of smooth objectives, it is known that such 
Coordinate Pivot Iterations (CPI) schemes converge to local minima (see for example 
[4j, p272). Further, in certain special cases of non-smooth objectives, it is also known 
that CPI schemes converge to local minima (see for example BZD- But in our case, the 
non-separable piecewise structure of the objective often causes our heuristic to halt at 
a point that is not a local minimum. Nevertheless, the global search heuristic is fruitful 
when utilized in a combined local-global search heuristic. This heuristic performs global 
searches with different initial points, each followed by a local search. We present nu¬ 
merical evidence, illustrating that it performs extremely well. Often finding the global 
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optimum in very few steps. 

The structure of the sequel is as follows. In Section [2] we present the model and 
basic properties. In Section [3] we focus on arrival departure dynamics, showing a piece- 
wise affine relationship between the arrival and departures times. We give an efficient 
algorithm for calculating the departure times given arrival times or vice-versa. This 
also solves the scheduling problem for the special case 7 = 0 . In Section [4] we charac¬ 
terise the constraints associated with quadratic programs which make up the piecewise 
quadratic cost. These are then used in the exhaustive search algorithm. We then 
present the local search algorithm and prove it always terminates at a local minimum 
(of the global objective). In Section Owe present our global search method based on 
CPI. We utilize the structure of the problem to obtain an efficient single coordinate 
search within the CPI. Then in Section 0 the local search and global searches are com¬ 
bined into a unified heuristic. We further illustrate the power of our heuristic through 
numerical examples. We conclude in Section O Some of the proofs are deferred to the 
appendix. 

Notation: We denote x Ay and x V y to be the minimum and maximum of x and 
y , respectively. We define any summation with initial index larger than the final index 
to equal zero (e.g. J^i= 2 a i = 0)- Vectors are taken as columns and are denoted in 
bold. 1 E R^ denotes a vector of l’s and e ; E R w denotes a vector of zeros in all but 
the Pth coordinate, which equals 1. The indicator function is denoted by 1 . 


2 Model 

Our model assumes that there is a fixed user set A f = {1,... ,N} where the service 
requirement of each user, £, is the same and is set to 1 without loss of generality 
(this can be accounted for by changing the units of [3 and a). Then the equations 
determining the relationship between the arrival times vector a = (aq,. .., aw)' and 
the departure times vector d = (d ±,..., djy)' are 

/ di 

v(q(t))dt , where q(t) = £i{‘ 6 [%.<%]}■ ( 2 - 1 ) 

» j&M 

Using the linear slowdown rate function, oa, the equations are represented as, 

1 = f (3 — a( t{t E [aj, dj]} — l)dt, i — (2.2) 

These N equations can be treated as equations for the unknowns d, given a or vice- 
versa. We assume N < /3/a + 1 so that it always holds that v(q(t )) > 0. 

The cost incurred by user i is, 

Ci(ai, d/) = (di - d*) 2 +7 (di - a t ), (2.3) 
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and the total cost function, which we seek to minimise, is 

c( a ) = J^Ci(ai,di(a)). (2.4) 

i£Af 

We assume (without loss of generality) that the ideal departure times, d* = (d*,..., d* N )' 
are ordered, i.e. df < ... < d* N . 

Remark For clarity of the exposition we choose the cost, (12. -I ll to be as simplistic as 
possible. Practical straightforward generalizations to the cost and to the associated 
algorithms and heuristics are discussed in the conclusion of the paper. These include 
other convex penalty functions, ideal arrival times and a potentially different penalty for 
early and late departures. Our algorithms, can all be adapted for such cost functions. 

We first have the following elementary lemmas: 

Lemma 2.1 Assume that the arrivals, a, are ordered: a\ < 02 < ... < a n, then the 
departures, d ; follow the same order: d\ < d 2 < • • • < d/v- 

Lemma 2.2 For any a there is a unique d and vice-versa. 

As a consequence of the assumed order of d* and of the above lemma we assert that an 
optimal schedule can only be attained with an ordered a whose individual coordinates 
lie in a compact interval, as shown in the following lemma. 

Lemma 2.3 An optimal arrival schedule satisfies a < cq < ... < on < a, where 

,* N _ „ N 

1 P-a(N-l)’ N /3 — a(N — 1) 

We may thus define the search region for the optimal schedule: 

1Z = {a G R a : a < ai < ... < Oat < a}, 

and take our scheduling problem to be min aGK c(a). 

No strict condition on the joint order of a* and di can be imposed except for the 
requirement that a* < di for any i (the sojourn time of all users is strictly positive). 
We are thus motivated to define the following for i 6 A f: 

ki := max {k E Af : a*, < di } = min [k E Af : ak+i > dfi\, (2.5) 

hi := min [h E Af : dh > aj} = max [h E Af : dh -i < a*}- (2.6) 

The variable /q specihes the interval [a^ja^+i) in which di resides. Similarly the 
variable h l specihes that a, lies in the interval (d^-i, dh J. Note that we dehne a o, d 0 := 
—oo and ri jV +i, djv+i := oo. The sequences ki and hi satisfy some basic properties: 
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(i) They are non-decreasing and are confined to the set A f. (ii) From the fact that 
di < di we have that i < k^. (iii) Since d is an ordered sequence and also a* < di we 
have hi < i. (iv) We have hi = 1 and k N = N. (v) Each sequence determines the 
other: 

ki = max {k E J\f : hk < i}, and h l = min {h G Af : k h > i}. 

Thus given either the sequence ki, i G Af or the sequence hi, i 6 Af or both, the 
ordering of the 2N tuple (ai,..., ajv, d\,..., d^) is fully specified as long as we require 
that ads and di s are ordered so as to be consistent with Lemmas 12.11 and [2731 
We denote the set of possible k = (ki, ..., k m)' by 

K := {k e Af N : k N = N, ki < k, \/i < j) . (2.7) 

Similarly, we denote the set of possible h = (hi,..., h N )' by Tt. We have that, 


\K\ = \n 


JV + l' 


This follows (for example) by observing that the elements of /C correspond uniquely 
to lattice paths in the N x N grid from bottom-left to top-right with up and right 
movements without crossing the diagonal. The number of such elements is the iV’th 
Catalan number, see for example p259 in [13]. 

The following example illustrates the dynamics of the model (without optimiza¬ 
tion) and shows the role of k, or alternatively h, in summarizing the piecewise affine 
dynamics. 


Example 2.1 Take (3 = 1/2, a = 1/6 and N = 3. This 3 user system exhibits rates 
that are either 1/2,1/3 or 1/6 depending on the number of users present. The free flow 
sojourn time is 1/(3 = 2 . Assume a\ = 0, a 2 = 1 and a 3 = 3. We now describe the 
dynamics of the system. See also Figure (M 

During the time interval [0,1), q(t) = 1 and the first user is being served at rate 
1/2. By time t = 1 the remaining service required by that user is 1/2. At time t — 1, 
the number of users in the system, q(t), grows to 2 and the rate of service to each user 
is reduced to 1/3. This means that without a further arrival causing further slowdown, 
user 1 is due to leave at time t = 2.5. Since 2.5 < a 3 , this is indeed the case. At 
t = 2.5, q(t ) changes from 2 to 1. By that time, the remaining service required by user 
2 is 1/2. Then during the time interval [2.5, 3) user 2 is served at rate 1/2 reducing the 
remaining service of that user to 1/4. At time t = 3, user 3 joins, increasing q(t) back 
to 2 and reducing the service rate again to 1/3. User 2 then leaves at time t = 3.75 
and as can be verified using the same types of simple calculations, user 3 finally leaves 
at time t = 5.25. 

Observe that for this example, the order of events is: 


a 3 < a 2 < di < a 3 < d 2 < d 3 . 


This then implies that for this schedule, 

k\ =2, k 2 = 3, k 3 = 3, and hi = 1, h 2 = 1, h 3 = 2. 
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Figure 2: An illustration of the dynamics of a three user example. The shaded gray 
areas show the remaining work for each individual user. Work is depleted at rate | 
when only one user is present and is depleted at the slower rate of | when two users 
are present. 


3 Arrival Departure Dynamics 

We now investigate the relationship between arrivals and departures, induced by the 
linear slowdown dynamics. 

Proposition 3.1 Equation (\2.2\i can be expressed as 

z—1 ki 

(f3 — a{ki — i))di — a dj — (j3 — a(i — hi))ai + a aj = 1, i G A7, (3.1) 

j=hi j=i+1 


or alternatively, 


D d — A a = 1, 


(3.2) 


with the matrices A 
' 

j3 — a{i — 





e 

and D e R w defined as follows: 






hi), 

i = 3, 

fi - a(ki - i), 

* — 3i 


i + 1 < j < ki, Dij \= < 

—a, 

K < j <i — 


o.w. 

o, 

o.w. 
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Proof We manipulate (12. 2 p to get, 




A f d 

'i 


(/3 + 

a)(dj 

- cii) - a 22 / 

j — 1 d a i 

l{i 6 [aj, dj]}dt 


09 + 

a)(dj 

N 

- Cli) - OL ^2 (di 

7 — 1 

A dj — diV a,j) + 




J — x 

2—1 


N 

(/3 + 

a)(dj 

cij) o ^ ] (dj 

A dj — aj V aj) + — a(dj — aj) 

— a (d. 



j=i 


i=»+i 



«— i 

N 


/3(di 

- a*) - 

- a (dj A dj - 

- a,i V Oj) + — a (dj A dj - 

- a* V aj) + . 



3 = 1 

j=i +1 



where in the third step we have used the fact that a* < dj for the term corresponding 
to j = i. We now use the fact that a and d are both ordered to get, 

i -1 N 

1 = f5(di — at) — a (dj — a*) + — a (dj — aj) + 

j= i j=i+ 1 

i— 1 N 

= /3(di — cii) — a dj — di A dj) — a (dj — aj A dj) 
t=i j'=»+i 

2—1 2—1 N 

= —fdai + (f3 — a(iV — i))dj — a dj + a A dj) + a (aj A dj). 

t=l j =1 j'=i+l 

Now the summations X^j=i( a * ^ dj) and XljLj+i( a i ^ can t> e broken up as follows: 

2 — 1 2 — 1 2 — 1 

y~~^(aj A dj) = l{dj < aj}dj + l{dj > ajja, 

j = 1 j=i j=i 

—1 2—1 /lj — 1 

^ ^ dj + ^ ^ Oj% ^ ^ dj H“ (z 
j= i j=hi j= i 

N N N 

(aj A di) = l{aj > dj}dj + l{aj < dj}aj 

jsi+l j=j+l j=i+l 

iV ki ki 

^ ^ dj ~t~ ^ ^ flj (W fcj)dj d" ^ ^ Uj. 
j=fci+l j=*+l j=j+l 
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Combining the above we obtain: 


i—1 hi — 1 ki 

1 = — (f3 — a(i — hi))ai + (j3 — a(ki — i))di — a ^ dj — d^j + a aj 

j =i j =i i=*+i 

i—1 ki 

— — ((3 — a(i — hi))ai + (j3 — a(ki — i))di — a ^ dj — aj j. 

j=hi j=i +1 

Rearranging we obtain (13.11) . §j 

The following observations are a consequence of Proposition 13.11 


1. Consider some user i arriving at time a* to an empty system, and departing at 
time di to leave an empty system. In this case there are no other users effecting 
his sojourn time or rate. For such a user k t = h, = i. In this case (13.11) implies 
that di = di + 1//3 as expected. 

2. The matrices A and D are lower and upper triangular, respectively, with a non¬ 
zero diagonal, and are therefore both non-singular. 

3. For the special cases i — 1 and i = N (using the fact hi = 1 and kjsr = N): 


k i 

(/3— a{k\— l))rfi —(3 ai+a a,j = 1, 

1=2 


N -1 

and /3dN~ct dj —(/3—a:(.ZV—/ ijv))ojv 

j=h N 


1 . 


I.e., 


1 + /3ai — a Ej =2 a j 
(3 — a (hi — 1) 


CLN 


Pd N - a Y,?=h N d j ~ 1 
(3 — a(N — h N ) 


The above structure suggests iterative algorithms for either determining a based 
on d or vice-versa. In both cases, k and h are found as bi-products. As an aid to 
describing these algorithms, define for i,k,h G A f and for a given a (respectively d), 
the functions di t k,h{- | a), a^hi • I d) : ~R N —* R as follows, 


d 


i,k,h 



&i,k,h 




1 + (p-a{i- h))di + a ( Y,]=h dj ~ Ej= i+ i a i) 
f3 — a(k — i ) 

(P ~ a ( k - i))d* - a( E 7=h d j ~ E?=i+i - 1 

/3 — a(i — h) 


Observe that in the evaluation of these functions, the arguments, d or a are only 
utilized for the coordinates indexed h ,..., i — 1 or i + 1,..., k respectively (if i — 1 or 
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respectively i = N these index lists are empty). Further observe that stated in terms 
of d(-) or a(-) and given k £ /C and h £ %, equations (]3.1|) can be represented as, 

d{ • • • 5 djy'j | . 

••, a N ) j, 

i e JV, 

or alternatively, 



CLi ( (^1? • • • 5 ^iv) | (^1? • 

• •, d N ) J, 

i £ A f. 


Given the above we have two (dual) algorithms for determining the network dynamics. 
Algorithm la finds the departure times based on arrival times. Algorithm lb finds the 
arrival times given the departure times. 

Proposition 3.2 Algorithm la finds the unique solution d to equations (12.2ft . given a. 
Similarly Algorithm lb finds a unique solution a to the equations, given d. Both algo¬ 
rithms require at most 2N steps in each of which (13. ip is evaluated. 


Algorithm la: Determination of network dynamics with given arrival times 
Input: a £ R N such that ay < 02 < ... < ajsr 
Output: d = (d\, ..., d/v), k = (Ay,..., kjy) and h = (hi, ..., hjv) 
init k = h = (1, 2, 3,..., N ) 
init d = 0 

for % — 1,..., N do 

set k = i V ki-i (taking k 0 := 1) 
compute dfik, hi, d\ a) 
while dfik, h, d | a) < ak +1 do 
increment k 
compute dfik, hi, d | a) 
end while 
set ki = k 

set di = dfik, hi, d | a) 

set h i+ 1 = max {h £ {1 ,..., i + 1 } : kh > i + l} 

end for 

return (d, k, h) 


3.1 Optimizing for Extreme Cases of 7 

As described in the introduction, optimizing ()2.4j) when 7 = 0 or 7 00 can be done 

efficiently. For the case 7 = 0 , all that is needed is to schedule arrivals so that each 
departure time, d* is exactly at d*. This achieves zero costs. Such a schedule is simply 
obtained by running Algorithm lb with input d = d*. This immediately leads to the 
following corollary of Proposition 13.21 
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Algorithm lb: Determination of network dynamics with given departure times 
Input: d G M. N such that di < d 2 < ... < 

Output: a = (ai,ajv), k = (k \,..., k^) and h = (hi, ..., hjy) 
init k = h = (1,2,3,..., N ) 
init d = 0 

for i = N,, 1 do 

set h — i A h l+ \ (taking hw+i N) 
compute a.i(ki, h, a | d) 
while a>i(ki, h, a | d) > d^-i do 
decrement h 
compute di(ki, h, a | d) 
end while 
set hi = h 

set a,i = a.i(ki, h, a | d) 

set ki -1 = min [k G {i — 1 ,..., N} : hk <i — l} 

end for 
return (a, k, h) 


Corollary 3.3 For the special case 7 = 0 there is an efficient polynomial time algo¬ 
rithm that finds the unique optimal schedule, a 0 , achieving c(a°) = 0. 


For the case of large 7 it is sensible to consider a classic schedule where users do not 
overlap: 



di<a i+ 1 , i = l,...,N-l. 


(3.3) 


This poses the problem as a classic single machine scheduling problem with due dates 
(see for example [3] or 12a). This implies that the total costs due to sojourn times is at 
the minimal possible value 7 N/(3 and the costs due to deviations from ideal departure 
times is, 

y^( a » + l//^ — d*i ) 2 ■ 

ieM 


For any finite 7 this does not necessarily minimize (12.411 . but as 7 —» 00 it is a sensible 
approximation. I.e. for large 7 the optimal schedule is approximated by the solution 
of the following convex quadratic program: 


min 

S.t. 


N 


Ytoi + 1//? - d*) 2 


i— 1 


ai-a i+1 <--, i — 1,... ,N — 1. 


(3.4) 


As for the case 7 = 0 , the above quadratic program can be efficiently solved using 
any standard convex quadratic programming method. Denote the optimizer by a°°. 
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3.2 A Linear Approximation 

Having the schedules a 0 and a°° for the cases 7 = 0 and 7 = 00 respectively, we are 
motivated to suggest a set of potential (initial) guesses for the optimal schedule for 
arbitrary 7. Let M > 1 be some integer specifying the number of initial guesses. Then 
the set of initial guesses lie on the segment interpolating a 0 and a°°: 

A ={ a ° jffj+ a “ i 1 ~ wfi) ■ m = 0 "-- Af — 1 }- (3 - 5) 

when M > 2 or equals {a 0 } if M = 1. We shall use the M points of A as initial 
guess points for the optimization heuristics that we present in the sequel. This is 
a sensible choice since every set of due dates d\,...,d* N exhibits some contour in 1Z, 
parametrized by 7, corresponding to the optimal schedules (for each 7). The end points 
of this contour are a 0 and a°° which we can efficiently find. Thus for a G [ 0 , 1 ], the 
points a 0 a + a°° (1 — a) constitute a linear approximation of this contour. In cases 
where the contour is almost not curved we have that the optimal value lies very near 
to the linear approximation. I 11 other cases, this is simply a set of initial guesses, yet 
possibly a sensible one. Note that the values of M do not need to be excessively large 
because initial points that are close are likely to yield the same local solutions. The 
numerical analysis of Section [6] reinforces this observation. 

4 Piecewise Quadratic Formulation 

Our key observation in this section is that the search region 1Z can be partitioned into 
polytopes indexed by k G /C, where over each such polytope, the objective is of a convex 
quadratic form. This yields \IC\ convex quadratic problems, each of which (individually) 
can be solved efficiently. An immediate exhaustive-search algorithm is then to solve 
all of the problems so as to find the minimising one. This yields a finite-time exact 
solution and is a sensible choice for small N (e.g. N < 15). But since, 

4 v 

1^1 ~ AT3/2 0F’ 

solving all convex problems is not a viable method for non-small N. We thus also spec¬ 
ify a local-search algorithm which searches elements of fC by moving across neighbouring 
polytopes until finding a local optimum. 

The following proposition is key: 

Proposition 4.1 The region 1Z can be partitioned into polytopes indexed by k G 1C, 
and denoted 

Pk := {a G TZ : a ki < [0 k a + r] k ]i < a ki+1 , i G W} , 

where @ k = D~ X A and rj k = ZW 1 1 with A and D based on k are specified by Proposi¬ 
tion EH Then for a G P k the objective function is convex and is given by, 

c k (a) = a'Q k a + b k a + b k , 
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with, 


Qk ~ © k ©k, 

b k = 2 (r 7 k -d*)'© k + 7 l'(© k -/), 

^k = (r/ k - d*)'(r/ k - d*) + rj k . 

Proof The results of Proposition 13.11 show that every k e /C specifies matrices D and 
A such that, d = D~ l A a + D~ l 1 = @ k a + r/ k . This holds with constant © k and 
77 k for all a and d for which k as defined in (12.51) is fixed. The polytope P k specifies 
this exactly by describing the set of arrival points for which the specific ordering of 
departures within arrivals is given by k. 

Since for all a e P k the affine relationship between a and d holds with the same 
© k and 77 k the cost, (j2.4j) . can be explicitly represented in terms of a: 

c(a) = - d*) 2 + 7 (di - a* ) 

ieAT 

= (d — d*) 7 (d — d*) + 7 1 / (d — a) 

= (a'0 k + ( 77 k ~~ d*)') (0 k a + 77 k — d*) + 7 l' (0 k a + rj k — a) 

= a'@k@ k a + (2(77 k - d*)'©k + 7 l'(0 k - /))a + (rj k - d*)'(r] k - d*) + ^l'rj k . 

This yields Q k , b k and the constant term, h k . Finally, since Q k is a Gram matrix, it 
is positive semi-definite. Hence the objective is convex. | 

4.1 Exhaustive Search 

We are now faced with a family of convex quadratic programs. For each k e /C, denote 
Ck(-) to be the cost associated with k then, 

QP( k) : minc k (a). (4.1) 

aGP k 

Note that while the constant term h k is not required for finding the solution of 
QP( k), it is needed for comparing the outcomes of the quadratic programs associated 
with different elements of JC. Indeed the most basic use of QP( k) is for an exhaus¬ 
tive search algorithm which finds the global optimal schedule in finite time. This is 
summarised in Algorithm 2. 

The virtue of Algorithm 2 is that it finds the optimal schedule in finite time. But 
this is done by solving an exponential (in N) number of convex QP(-) problems, so for 
non-small N it is not a sensible algorithm. Hence we now introduce a search heuristic. 

4.2 Neighbour Search 

In this section we introduce a heuristic search aimed at finding a local minimum by 
searching on neighbouring regions. The search procedure solves the QP (14.11) over 
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Algorithm 2 : Exhaustive search for global optimum 
Input: Model parameters only (N,a,{3, d* and 7 ) 
Output: a* (global optimum) 
init m* = cxd 

for k G /C do 

solve QP( k) with optimise!' a and optimum m 

if m < m* then 

set a* = a 
set m* = m 

end if 
end for 

return (a*, m*) 


neighbouring elements of /C by changing a single coordinate of k at a time. We prove 
that this procedure converges to a local minimum; yet this may possibly take an ex¬ 
ponential number of steps in the worst case. 

Given a solution a of QP(k ) we define the following two sets of indices: 

Xi(a, k) := {j G N : [0 k a + r] k ]j = a kj+1 }, 

X 2 (a, k) := {j G Af : a kj = [0 k a + rj^}. 

Noting that d* = [0 k a + r] k ]i, and recalling that k t is index of the maximal arrival 
time that is less than or equal to d* we have that if i G X\ (a, k) then the optimal 
solution of QP{k) exhibits d, : = a ki+ 1 as an active constraint. Hence a neighbouring 
region to the constraint set P k is P k (i> where k ( d = k on all coordinates except for 
i where it is equal to hi + 1. Similarly if i G X 2 (a, k) then a ki = d* as an active 
constraint. In this case, k ( d is set to equal k on all co-ordinates except for i where it 
is set to equal ki — 1. Thus for every element of Xi(a, k) and X 2 (a, k) we have a well 
defined neighbouring region. Defining now the sets of neighbouring regions to P k by 

JC t {h{ a, k)):={k« : iGl^a, k)}, £ = 1 , 2 , 

we have the following local search algorithm: 
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Algorithm 3 : Neighbour search for local optimum (local search) 

Input: k 

Output: a* and m* 

solve QP( k) with optimiser a and optimum m 

init m* = m 

init a* = a 

for i E /Ci(Xi(a, k)) do 

solve QP( k ( ^) with optimiser a and optimum m 
if m < m* then restart algorithm with k — k !l) 

end for 

for i G /C 2 (X 2 (a, k)) do 

solve QP( k^) with optimiser a and optimum m 
if m < rtf' then restart algorithm with k = k (n 

end for 

return (a*, rrf) 


Proposition 4.2 Algorithm 3 converges to a local minimum for any initial vector k. 

Proof Every step of the algorithm can only improve the objective function, since 
m < m* is the condition for the change of k, hence the algorithm cannot go back to 
a region which it has already visited. Furthermore, there is a finite number of regions 
which means the algorithm terminates in a finite number of steps. If for some a which is 
the solution of QP( k) there are no improvements in any of the neighbouring regions the 
algorithm stops at a local minimum. This can be either due to no active constraints to 
QP( k) (an interior point) or due to the fact that the neighbouring quadratic programs 
do not improve on the solution of QP( k). | 

5 Global Search Over Single Coordinates 

In this section we put forward Algorithms 4 and 5 that together form a coordinate 
pivot iteration procedure. We first describe how the dynamics presented in Sections [2] 
and [4] can be used to find a global minimum with respect to a single coordinate r 6 J\f 
(user) when all other coordinates are fixed. We call this procedure a global search over 
a single coordinate r. 

The computational complexity of such a procedure is shown to be at most 0(N 5 ). 
We then utilise this procedure to define a coordinate pivot iteration algorithm, that 
performs optimization cycles on all of the coordinates until no improvement can be 
made. 

To understand the main idea consider Figure [3] (a). This figure corresponds to an 
example with N = 4, a = 1.5 and /3 — 5. Here the arrival times a 2 , 03 , <24 are fixed 
at (0.05,0.15,0.45) and the arrival time of user r = 1 (denoted also x) is allowed to 
vary. The (horizontal) blue dotted lines denote the fixed arrival times a 2 , 03 , 04 . The 
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thin blue curves correspond to the departure times d 2 , d 3 , (I 4 . The thick green dotted 
and solid curves correspond to the arrival and departure time of user 1 respectively. 
When x is small enough or large enough, it is seen that user 1 does not affect the other 
users. But otherwise, user 1 interacts with the other users and potentially modifies 
their departure times. 



d\ 


d\ 


d 


4 

4 

3 

2 


0i 3 

d2 


x = a 1 



Figure 3: (a) Arrival (horizontal dotted) and departure (horizontal solid) profiles ob¬ 
tained by changing the arrival time of user 1 . (b) Cost function obtained by changing 
the arrival time of user 1. Break points are marked in both (a) and (b) by vertical lines 
as follows: Solid black lines mark Type la points (note there are exactly N — 1 = 3 
such breakpoints). Dotted black lines mark Type lb breakpoints (note that there are 
exactly N — 1 = 3 such breakpoints as well). Type 2a breakpoints are marked by 
dashed red lines and Type 2b breakpoints are marked by brown dashed-dotted lines. 


As is further evident from Figure |3] (a), the dynamics of the departure times are 
piecewise affine with breakpoints as marked by the vertical lines in the figure. In 
between these lines, the effect of changing x on other quantities is affine. In between 
these breakpoints, the objective function is piecewise convex (quadratic). This property 
is illustrated in Figure [3] (b) where the objective is plotted as a function of x. This 
property allows us to optimise globally over a single coordinate, utilizing the problem 
structure. The desired departure times used for the cost function in (b) were d* = 0.5 
for i = 1,..., 4. 
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The global search over a single coordinate works by varying x from a to a and in the 
process searches for the one-coordinate optimum. This is done with a finite number of 
steps because of the piecewise-affine dynamics. Our algorithm incrementally computes 
the piecewise-affine coefficients within these steps. We call each step a “breakpoint”. 
The following types of breakpoints may occur: 

Type la: The arrival of r overtakes the next arrival of any i 
(solid black line). 

Type lb: The departure of any i is overtaken by the arrival of r 
(dotted black line). 

Type 2 a: The departure of any i overtakes any arrival 
(dashed red line). 

Type 2 b: The departure of any i is overtaken by an arrival of j ^ r 
(brown dashed-dotted line). 

Observe that in varying x, breakpoints of type la and lb occur exactly N — 1 times 
each. Less trivially, we have a bound on the number of type 2a and 2b breakpoints: 

Proposition 5.1 In executing the global search over a single coordinate r, the total 
number of breakpoints is 0(N 3 ). 

Before presenting the proof, we present the details of the piecewise-affine dynamics and 
the details of the global search over a single coordinate r algorithm. 

5.1 Algorithm Details 

In carrying out the global search over a single coordinate r, we remove the restriction 
that arrival times are ordered. That is, the search region is extended from IZ to a set 
not requiring such order IZ := [a, a} N . This allows us to carry out a full search for the 
optimum with respect to a single user r without the restriction a r G [a r _i, a r+ \}. This 
broader search potentially enables bigger gains in the objective when integrating the 
algorithm within a search heuristic. Further, any point a G IZ can be mapped into a 
unique point 0( a) G IZ where O(-) is an ordering operator. By Lemma [2.31 we have 
that c((9(a)) < c(a). 

Take a G IZ as an initial arrival vector and suppose that we are optimising over 
user r. Let x G [a, a] be the immediate search value of a r (keeping the other arrival 
times fixed). For any such x we define a corresponding permutation 7 r(a, x) indicating 
the current order of arrivals, as well as the ordered arrival vector 

a(a, x) . (D (o^i ( x ) j • • • j ^■n r {x)—l i ®7r r (:r)+l j • • • j ®7Tjv(x)) ■ (5.1) 


18 




This vector can serve as input to Algorithm la yielding a corresponding d(a, x), k(a, x) 
and h(a, x). Furthermore, using (13.11) we have the local piecewise-affine relationship, 


G?j(o, x) X 9i 1^77(3, x),k(a, x) "F 7ji |r,7r(a, x),k(a, x) j ^ F A/", X G [o, Cl]. 


That is, the coefficients of the departures between breakpoints depend on the permu¬ 
tation of the users as well as on the current order of their arrivals and departures. For 
brevity we omit the dependencies on x. a, 7r and k. Manipulating (13.ip we obtain, 


(0i, Vi) = 


0 , 


l+a i {P-a(i-hi))-a(j 2 l jLi + i a j-T,}=h i 




9 


0 +-E) 




\ P-a{ki-i) 


l+°* (P— a (i—hi))—oiajl{j^n r }-E 
’ /3-a(ki-i) 

i-«(E*l <+1 °,-E 

f3—a(ki—i) 




l < 

7T r 5 ‘ 

fcj < 7T, 

i < 

7T r , 1 

ki > 7T, 

i = 

7T f 7 



aE'"l On 
/3-a(ki—i) ’ 


l+ai(p-a(i-hi))-a(j2 k jL i+ i a j~Yl)Jh % 


f}—a{ki—i) 


l > 7l r 


(5.2) 

On every interval, the departure times di are all affine and continuous w.r.t x with 
the above coefficients, until a breakpoint (of type la, lb, 2a or 2b) occurs. Comput¬ 
ing the time of the next breakpoint is easily done by considering the piecewise affine 
dynamics. Potential breakpoints of types la and lb are to occur at times t where 
x +1 — a nr+ 1 and t dj + di = a r + t, respectively. Potential breakpoints of types 2a and 
2 b involving user i are to occur at times tdi + di = a^+i and tdi + di = a.k, respectively. 
Observing now that type 2a breakpoints may occur only when di > 0 and type 2b 
breakpoints may occur only when di < 0 we have that the next breakpoint occurs at, 


r = min{« 0 , t u .. ,,t N ,t N+1 }, 


(5.3) 


where t Q = 
1 < i < N: 


a n , r+ \ — x (type la breakpoints), t N+ 1 = a — x (termination) and for 


U 


a k .-9iX-r]i 

el 

a k —0iX—r)i 
Oi-1 

a k.i+i 9{X ik 

el 

00 

< 


, 0 i < 0 , hi ^ r, 

, 0 i < 0 , hi = r, 

, di > 0, 

, 0 i = 0. 


Considering the time interval until the next breakpoints, [x, r] we have that the total 
cost as a function of the arrival time x G [x, r] of user r is 


c(x; 7r) := ^ ((^.x + r) v . - d*f + 7 {d n .x + rj Vj - a nj )) , 
j€ M 
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with derivative <9c(x; 7 r) = JT<=a r 0^ j (2(r) nj - d*) + 7 ) + 2 x and with the root 

Xo > x, solving dc(x 0 ; 7r) = 0 (and often not lying within the interval [x, r]): 

_ ~ {‘ZjVirj ~ dj) + 7 ) 

X ° _ 2 V 0 2 

z 2^ jeJV % 

Note that it is crucial to keep track of 7r at every step in order to associate the correct 
ideal departure time to every user. In iterating over intervals we search for the minimal 
c(x; 7r) (denoted m*) as follows: If <9c(x; 7r) > 0 for all x e [x,r], then we continue to 
the next interval. Otherwise, if x 0 — x < r and m* > c(x 0 ; 7r), then set m* = c(x 0 ; 7r) 
and (a*) = xo, and if Xo — x > r and m* > c(x + r; 7r), then set m* = c(x + r; 7 r) and 

(a*) r = x + r. 

In this way, x updates over intervals, of the form [x, r\. Prior to moving to the 
next interval we need to update the permutation variables 7 r, k, and h. Denote the 
minimizing set of (15.3ft by T := argmin{f 0 , h ,..., tjsr} and sequentially for every i G T: 

• If i = 0, then we update 7 r by changing the order between user r and the next 
user j : 77 = 77 + 1, i.e. set 77 = 77 + 1 and 7 Tj — iTj — 1. In this case, there is 
no change in k or h. (Type la breakpoints). 

• If i 6 {1,..., N}, then the order 7 r does not change, but we update k and h: If 
9i < 0, then update h k . = h ki + 1, followed by ki = ki~ 1. If 0* > 0, then update 
ki = ki + 1, followed by h ki = h ki — 1. (All other types of breakpoints). 

• If i — N + 1, then the iteration is complete and no changes are required. 

Remark For any convex and differentiable cost functions, the first order condition 
yielding Xo can be solved. For some elaborate functions this may also require a numer¬ 
ical procedure. If the late and early cost functions are not strictly convex (for example 
affine), then computing xo can be skipped. If the cost function is piecewise affine, then 
only the sign of dc needs to be computed, and if it is negative check if the next point 
x + r is a new minimum point or not. 
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Algorithm 4: Global search over a single coordinate 

Input: a G U, r e AT, 7T 
Output: a* and m* 
init x = a r = a 
init a = 0( a) 
run Alg.la(a) — > (d, k, h) 
init a* = a 
init m* = c(x; tv) 
set 7 T r = 1 
for i < r do 

Set 7 Ti = 7Tj + 1 
end for 
while x < a do 
set a = 0 (a) 

compute: 6, 77 , r, T, and dc(x; tv) 
if dc(x; 7r) < 0 then 

compute xq and c(x 0 ; 7v) 
if xo < x + t then 

if c(x 0 ; tv ) < m* then 
set a* = x 0 
set m* = c(x 0 ; tv) 
end if 

else if c(x + r; tv) < m* then 
set a* = x + r 
set m* = c(x + r; tv) 

end if 
end if 

set x = x + t 

for i G T do 
if ? = 0 then 

set 7 Tj = TVj — 1 where j satisfies — Tv r + 1 

Set 7T r = 7T r + 1 

end if 

if i G {1,..., N} then 
if 9i < 0 then 

set h ki = h k . + 1 and hi = hi — 1 

else if 9i > 0 then 

set hi = hi + 1 and h ki = h ki - 1 

end if 
end if 
end for 
end while 
return (a*, m*) 
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5.2 Computational Complexity 

In the following series of lemmata we analyse the complexity of Algorithm 4. In 
particular, we prove Proposition 15. 11 establishing bounds for the number of breakpoints 
of each type. Throughout the analysis we continue denoting the coordinate being 
optimised by r and the respective value by x = a r . Keep in mind that 7p, 9 t , ki , and 
hi are functions of x and the initial unordered vector a for every i G A f. We treat a as 
the ordered vector (15.ip as before. 

Lemma 5.2 For any i G Af such that * ^ r, the coefficient 9i < 0 and as a consequence 
dj(x) is monotone non-increasing for every x > ai. 

Lemma 5.3 For any permutation n at the start of the global search on r, the coefficient 
9i of any i G n changes sign from strictly positive to strictly negative or vice versa at 
most i — 1 times during the search. 

We now prove proposition 15.11 

Proof For any 2 < i < N in the original permutation 7r, the type 2a and 2b 
breakpoints occur at most N — i times for every change of sign. This is because their 
departure time can only cross arrival times of later arrivals. According to Lemma 15.31 
the number of sign changes for any 2 < i < N is at most i — 1. Thus, the total number 
of breakpoints of type (2a or 2b) is at most 



i=2 


Thus, adding up all types of breakpoints, we get that the search domain [a, a] is broken 
up to at most (| N 3 — N 2 + |N — 2) intervals. | 

Furthermore, we have the following bound for the complexity of Algorithm 4. 

Corollary 5.4 The computation complexity of Algorithm 4 is at most 0(N 5 ). 

Proof In every interval step of a global search on a single coordinate there is a need 
to compute the coefficient vectors r] and 6 . This is equivalent to calculating the 
departure times recursively using Algorithm la. In Proposition 13.21 it was shown that 
the recursion requires at most 2N steps. On top of this, in every one of these steps the 
actual computation requires summation of up to N variables. Now since the number 
of breakpoints intervals is bounded by 0(N 3 ) we conclude the result. J 


22 











5.3 Coordinate Pivot Iteration Optimization 

In this subsection we illustrate how Algorithm 4 can be applied to carry out standard 
Coordinate Pivot Iteration (CPI), see [4], p272. In every iteration of the CPI algorithm, 
the total cost function is minimized with respect to the arrival time of one user, when 
all other arrival times are fixed. This is then repeated for all users; we call the iteration 
over all N users a CPI cycle. The CPI algorithm stops when the total improvement in 
a cycle is smaller than some specified tolerance parameter, e > 0. Note that in non¬ 
smooth CPI (such as our case), CPI often stops when the total improvement is in-fact 
exactly 0. That is, e is often not a significant parameter. A further comment is that 
our CPI algorithm utilizes Algorithm 4 searching over the broader space, 1Z. We can 
thus improve the objective (see Lemma I2.3|l by incorporating the ordering operator, 
O , at the end of each CPI cycle. 

We add the following notations for the optimization procedure: Let n = 0,1,... 
be the cycle number, d n ' the total cost at end of cycle n, m* the global minimal total 
cost, and a* the global optimal arrival vector. 


Algorithm 5: Coordinate pivot iteration (global search) 

Input: a^ and e 
Output: a* and m* 

init n = 0 
init A = e + 1 
init a* = a^ 0 ) 
init d°) = c(a*) 
while A > e do 
set n — n + 1 
set a = a* 
for r 6 Ado 

run Alg. 4(r, a) —> a, 
end for 
set a* = 0( a) 
set c= c(a*) 
set A = c (n-1) — 
end while 
set m* = 
return (a*, m*) 


Hinging upon the results of the previous section, we have: 

Corollary 5.5 The computation complexity of a single CPI cycle, i.e. conducting a 
line search on all coordinates, is at most 0(N 6 ). 
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Proof In Proposition 15.41 we established that for a single coordinate the complexity is 
at most 0(N 5 ). It is therefore immediate that the complexity of running the algorithm 
for every coordinate is at most 0(1V 6 ). | 

Note that while we have a polynomial time CPI algorithm, there is no guarantee 
that it converges to a local minimum since the objective function is not smooth. In fact, 
numerical experimentation suggests that this is typically the case when the number of 
users is not very small, i.e., N > 4. Nevertheless, experimentation has shown that CPI 
algorithm generally outputs an arrival vector that lies in the vicinity of the optimum. 
This motivates combining it with the neighbour search, Algorithm 3 as discussed in 
the next section. 


6 A Combined Heuristic and Numerical Results 

We now utilise the problem structure and aforementioned algorithms to produce a 
combined heuristic. We use A as in (13. 5 p for initial points. For each of these M 
initial points we run a CPI (global) search followed by neighbour (local) search. The 
core principal is to use the CPI method in order to find a “good” initial polytope, or 
equivalently an arrival-departure permutation, and then to seek a local minimum using 
the neighbour search. 


Algorithm 6: Combined global and local search heuristic 
Input: Model parameters only (N,a,(3, d* and 7) 

Output: a* (local optimum) 
init m* = 00 
for aeido 

run Alg.5( a, ■ • •) —t (a,m) 
set k = k(a) 

run Alg.3(k, ...)—>■( a, m) 
if rh < 771 * then 

set a* = a and m* = rh 

end if 
end for 

return (a*, m*) 


We tested the combined heuristic Algorithm 6 on a variety of problem instances 
and it appears to perform very well both in terms of running time and in finding what 
we believe is a global optimum. Here we illustrate these results for one such problem 
instance. We take (3=1 and a = 0.8 /N (in this case the maximal slowdown is of the 
order of 80% independently of N ). We set d* as the N quantiles of a normal distribution 
with mean 0 and standard deviation 1/2. That is, there is an ideal departure profile 
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centred around 0. It is expected that when using optimal schedules, more congestion 
will occur as N increases and/or 7 decreases. 

Figure Q] illustrates the dynamics of the obtained schedules as generated by the 
heuristic (using M = 3 and e = 0.001). In these plots arrival times of individual users 
are plotted on the top axis, marked by blue dots, shifted to the right by the free flow 
time (1//3 = 1). Departure times are plotted on the bottom axis. Users that do not 
experience any delay are then represented by lines that are exactly vertical. Further, 
the more slanted the line, the more slowdown that the user experiences. The ideal 
departure times are marked by green stars. Hence ideally the stars are to align with 
the red dots. This occurs exactly when 7 = 0, and approximately occurs for small 7 , 
for instance 7 = 0.1 as in (a) and (d). Then as 7 is increased, the optimal schedule is 
such that there is hardly any delay (almost perfectly vertical lines), but in this case, 
users experience major deviations between departure times and the ideal values. 

For N = 15, as presented in (a)-(c), we were indeed able to verify optimality using 
the exhaustive search Algorithm 2. For N = 50, as presented in (d)-(f) we are not 
able to use the exhaustive search algorithm in any reasonable time. Nevertheless, in 
this case, in addition to seeing qualitatively sensible results, experimentation showed 
that increasing M does not modify the results. Hence we believe that the obtained 
schedules are also optimal. 

For N < 15, we were not able to find a case where the heuristic did not find the 
optimal schedule. This was tested on a wide range of parameter values by varying a 
and 7 and randomly generating multiple due date vectors. Further for large N (up to 
500) we see insensitivity with respect to M (the number of initial points) as well as 
to other randomized initial points. This result was also robust to changes in all of the 
parameter values (a , /3, 7 , and d*). This leads us to believe that our heuristic performs 
very well. 

Results, comparing running times are reported in Table [T| where we consider the 
algorithm with a single initial point a 0 (M = 1 ), and compare it to the exhaustive 
search given by Algorithm 2. For this table, we use the same problem data as described 
above, but scale the standard deviation by N to be 0.041V. For N < 15 the combined 
heuristic converged to the global optimum as verified by the exhaustive search with 
a negligible number of computations. For example, for N = 15 the heuristic method 
made ~ 737 core computations, i.e. solving a QP for a single CPI interval or NS 
polytope, in 3.28 seconds, while the exhaustive search had to solve ~ 10' quadratic 
programs and required about 11 hours Clearly, for larger N it is not feasible to run 
the exhaustive search while the combined heuristic is still very quick, as seen for up to 
N = 50 in Table [Q 

2 These computation times are using an AMD computer with 4 Phenorn II 955 3 . 2 GHz processors, 
with our algorithms implemented in version 3.1.2 of the R software. 
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Cti = • 


0 1 2 
(a) IV = 15, 7 = 0.1 


di = • 

d* = * 



(b) IV = 15, 7=1 



(d) N = 50, 7 = 0.1 




Figure 4: Optimal arrival-departure diagram for a = 0.8/IV, /3 = 1, and d* the N 
quantiles of a normal distribution with mean 0 and standard deviation 1 / 2 . 
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Table 1: Running time in seconds and computational steps of the combined heuristic 
(Algorithm 6 with M — 1) and the exhaustive search (Algorithm 2). 


N 

3 

5 

10 

11 

12 

14 

15 

20 

30 

50 

Combined heuristic 

CPI cycles 

3 

2 

4 

5 

4 

3 

3 

3 

3 

4 

Total breakpoints 

24 

65 

306 

382 

441 

642 

727 

1,383 

3,260 

8,636 

NS QPs solved 

2 

2 

9 

6 

8 

14 

10 

29 

34 

29 

Running time (sec.) 

0.05 

0.15 

1.17 

1.67 

1.99 

3.09 

3.28 

6.38 

18.31 

85.33 

Global opt. 

Yes 

Yes 

Yes 

Yes 

Yes 

Yes 

Yes 

NA 

NA 

NA 

Exhaustive search 

| /C | QPs solved 

5 

42 

16,796 

58,786 

208,012 

2.6 ■ 10 6 

9.7 ■ 10 6 

6.6 ■ 10 9 

3.8 ■ 10 15 

2 • 10 27 

Running time (sec.) 

0.00 

0.05 

25.25 

162.41 

509 

8,206 

39,454 

NA 

NA 

NA 


To further investigate our combined heuristic, in Figure [5] we illustrate the number 
CPI cycles and breakpoints, along with the respective number of quadratic programs 
solved by the neighbour search, until convergence of Algorithm 6 . The problem data 
was scaled as in the previous example. For every N the initial points given by A with 
M = 5 distinct initial points. The figure displays the minimum and maximum values 
out of the 5 initial points. Note that for every N the algorithm converged to the same 
local minimum for all initial points in A. 

We can see that the number of required CPI cycles was small and stabilized on 
2 regardless of the number of users. However, we should take into account that the 
number of coordinate iterations in every cycle is N, and that the complexity of each 
iteration also grows with N. Specifically, Proposition 15.41 shows that the number of 
breakpoints for every coordinate in the CPI is at most IV 3 , but in the example we see 
the growth is in effect linear (~ 3 N). Furthermore, the number of required quadratic 
programs solved in the neighbour search also grows linearly (~ |IV). This hints that 
the CPI does indeed find a point that is very “close” to a local minimum. The widening 
gap between the minimum and maximum number of NS iterations suggests that some 
of the initial points are better than others, and thus it is worthwhile trying several 
of them. The last point is important when solving for even larger values of N as the 
algorithm becomes more sensitive to “bad” initial points and may require setting a 
maximum number of iterations parameter for every initial point. Roughly, when 7 and 
a are both small, starting closer to a 0 is better and when they are both large, starting 
closer to a°° is better. However, for most combinations of parameters there seems to 
be no a-priori indication of what is a “good” starting point. Thus it is still beneficial 
to do the full search on A. Again we stress that the behaviour displayed in Figure [5] 
was robust with respect to changes in the model parameters. 
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(c) NS iterations 


Figure 5: Number of iterations in each component of Algorithm 6 as a function of N. 


7 Conclusion and Outlook 

We presented a model for a discrete-user deterministic processor sharing system, and 
addressed the problem of scheduling arrivals to such a system with the goal of minimiz¬ 
ing congestion and tardiness costs. A full characterisation of the congestion dynamics 
and an efficient method for computing them was provided. It was further shown that 
the optimal arrival schedule can be computed in a finite, but exponentially large, num¬ 
ber of steps. Several heuristics were therefore developed with the goal of an efficient 
computation of the optimal schedule. A combined global and local search heuristic 
was presented and numerically analysed. This method was shown to be efficient in 
numerical examples for a large population of users. 

The essential parts of our analysis and results applies for a much more general cost 
formulation, as we shall next detail. Given that user i enters the system at time a* and 
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leaves at time d* > a*, a plausible cost incurred by the user is the following: 

Ci{ai , di) = g- 1] ({di - d* ) + ) + gf ] ((d* - di) + ) (7.1) 

+ g? ] ((>* - a *) + ) + gf (« - ai ) + ) 

+ g? ] (di - a*), 

where (x) + := max(i, 0), and g!f\-),j G {1,... 5}, i G Af are some convex functions. 

The hrst and third terms of (17.11) capture the penalty for being late to the ideal 
departure and arrival times d* and a*, respectively. The second and fourth terms are 
the user’s cost for arriving and departing early. The fifth term is the user’s cost for 
travel/usage of the system. Our algorithm and results in this paper hold with slight 
technical modifications for arbitrary convex For purpose of exposition, we 

focused on, g^ 1 \x) = gf\x) = x 2 , gf\x) = gfHx) = 0 and gf\x) = 71. If adapted 
to the more general formulation, The exhaustive and neighbour search algorithms of 
Section |4] will generally require solving a constrained convex program, instead of convex 
quadratic, for every region. If g^\x) 7 ^ gf\x) and/or gf\x) 7 ^ gf 1 \x) 1 namely there 
are different penalties for arriving/departing later and early, then the CPI algorithm of 
Section 1531 will require some refinement of the definition of the piecewise segments. The 
complexity will not change as for every single coordinate there will be an addition of 
at most three segments, corresponding for these new points of discontinuity. Moreover, 
the root of the hrst order condition in every continuous segment will be given by the 
general form of the functions, instead of the quadratic root. 

An interesting generalization is considering a system with users who have hetero¬ 
geneous service demand. If this is the case then the order of departures is no longer 
identical to the order of arrivals. This means that the characterisation of Proposition 
13. 1 1 is no longer valid. 

A natural complementary model to this work is considering a decentralized decision 
framework in which the users choose their own arrival time. Namely, a non cooperative 
game with the individual arrival times are the actions of the players. This game is 
formulated and analysed in [ 20 ]. 

Finally, there is the challenge of characterising the computational complexity of our 
scheduling problem. We believe that finding the optimal k* G /C is an NP-complete 
problem but we still do not have a proof for this. Our belief is motivated (but not 
supported) by the fact that there are a number of related optimization problems which 
are known to be NP hard. Our problem is equivalent to a special case of one of them, 
namely non-linear integer programming. 

As we have shown, our goal is to minimize a non-convex piecewise quadratic ob¬ 
jective function, subject to piecewise linear constraints. It is known that non-convex 
quadratic programs and non-convex piecewise linear optimization are both NP hard 
(see [12] and |15j). In [23] it is shown that piecewise linear optimization problems can 
be modelled as linear mixed integer programs, where the definition of a piecewise linear 
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program relies on different coefficients for different polytopes, in a similar manner to 
our piecewise quadratic formulation in Section SJ It may be possible to apply similar 
methods with modifications for the piecewise convex instead of linear objective. How¬ 
ever, there is a more natural construction for our case. Let a(k) be the solution to 
QP(k), i.e., the solution to the local convex QP of a polytope k e 1C. But this can 
also be viewed as a function of the integer vector k which we can compute in polyno¬ 
mial time. Hence, solving our problem in polynomial time is equivalent to solving the 
non-linear integer program: 

min a(k)'Q k a(k) + b k a(k) + 6 k . 

ke/c 

Recall that fC = {k e J\f N : Py = N, k, < k 3 Vi < j } defines a set of linear con¬ 
straints on the integer decision variables. Clearly the objective is not linear with respect 
to k, as a(k) itself is already not necessarily linear. Such problems are known to be 
NP hard. See for example, [7\ and H3- Although we could not find a straightforward 
reduction of the problem to a known NP hard problem, we have shown that our prob¬ 
lem can be formulated as an (rather cumbersome) instance of a polynomial integer 
program, and have no reason to believe that the specific model comes with significant 
simplification of the general form. 

As a closing note we mention that it is generally of interest to compare our heuristics 
to potential integer programming methods. One may either discretize time and solve 
integer programs, or alternatively seek related integer programming formulations. It 
remains an open problem to compare our heuristics to such potential methods both in 
terms of accuracy and computation time. 
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A Proofs 


Proof of Lemma 12.11 Consider two arrivals a, < a ? . During the time interval [a,, a,-], 
user i has received some service, 



v{q(t))dt, 


while user j has not. Then during the time interval [ aj , d* A dj} both users receive the 
same service, J^ Adj v(q(t))dt. Then if d t > dj we have that J^ Adj v(q(t))dt = 1, which 
in turn would imply that, 



v[q(t))dt + 


cdiAdj 


>(q(t))dt 


**di 


v(q{t))dt > 1, 
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a contradiction. Hence di < dj. | 


Proof of Lemma 12.21 Without loss of generality assume a± < ... < a at and hence by 
the previous lemma, d is ordered. Assume now that there exists a d / d and define 
i = min{z : di ^ di}. Without loss of generality, assume that di < di. Using (12. ip it 
holds that, 

pdi pdi pdi 

/ v(q(t)) dt — 1 — / v(q(t)) dt+ v(q(t)) dt. 

J CLj j J Clj *J di 

Now since for all t < d t it holds that q(t) = q(t), then, 



v(q(t)) dt = 0. 


A contradiction. 


Now there exists a full symmetry between a and d, hence going in the opposite 
direction (for every d there exists a unique a) follows a similar argument to the above. 

I 


Proof of Lemma 12.31 We first argue that an optimal arrival must be ordered («i < 

... < a^v) by means of an interchange argument. Assume this is not the case, i.e. a is 
an optimal arrival schedule such that a* > dj for some i < j (such that d* < d*). If 
we switch between the arrival times of users i and j: di = dj and dj = di, while not 
changing any other arrival time, then because all users have the same service demand 
the departure times of all other users do not change. Consequently, the departure times 
are also switched: di = dj and dj = di. Therefore, the only change in the total cost 
function is the change in the cost incurred by i and j themselves. The change in the 
cost incurred by user i is given by (12.31) : 


Ci(di, di) - Ci(di, di) = (di - d *) 2 + 7 (d t - di) - (di - d *) 2 - 7 (di - di) 

= (dj — d *) 2 + 7 (dj — dj) — (di — d *) 2 — 7 (di — a*) ’ 

and for user j: 

Cj(dj, dj) — Cj(dj, dj) = (dj — d *) 2 + 7 (dj — dj) — (dj — d *) 2 — 7 (dj — dj) 
= (di - d *) 2 + 7 (di - di) - (dj - d *) 2 - 7 (dj - dj) 


Summing (1A.1D and (1A.2I) we obtain that the total change in cost is 


(A.l) 


(A.2) 


2 (d*-d*)(d i -dj). 

From Lemma \2. II we know that if di > dj then di > dj, and that by definition d* > d*, 
hence the change in the total cost function is negative which contradicts the assumption 
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that the schedule is optimal. In conclusion, any unordered schedule can be improved 
by a simple interchange of a pair of unordered coordinates, and therefore an optimal 
schedule must be ordered. 

The slowest service rate occurs when all N users are present in the system, and 
therefore the longest possible sojourn time is ^_ q ^ 1 jv _ 1) . The total time required to clear 

all users from the system is then coarsely upper bounded by . A schedule such 

that ai < a is clearly not optimal, since a trivial improvement can always be achieved 
by setting ci\ — a and shifting to the right the arrival times of any user that overlap 
due to the change in a\. We are guaranteed this is possible by the fact that all users 
can arrive and leave the system in the interval [a, d j], without any overlaps. Clearly, 
the deviation from ideal times can only decrease when making this change, while the 
sojourn times remain unchanged. The coarse upper bound a holds for the same reasons. 

I 

Proof of Proposition 13.21 The proof is for Algorithm la. The argument for Algo¬ 
rithm lb follows the same arguments. For every user i, iterating on all possible values 
of ki E J\f ensures that every possible departure interval [ak,ak+ 1 ) is checked. In a 
sense, this is an exhaustive search on all solutions that satisfy the dynamics given by 
Proposition 13. II Therefore, the algorithm will always converge to the unique solution. 

Given a vector of arrivals a E RA, for every i E A/", the departure time di occurs 
in one of the above defined partitions [ak,a,k+ 1 ), k E A f. The total number of steps 
will include the number of “correct” computations, that is for every i and ki = k the 
resulting di will indeed be in the interval [a^, a^+i). In total there will be exactly N 
correct computations. However, there will also be steps which will turn out to be false: 
for a given k the departure time di will not be in the interval [ak,ak+ 1 ). If kj = k for 
some j, then for every % > j: ki > k. Therefore, if for some i and k > i the computation 
will yield di [a*,, a^+i), then this interval will not be attempted by any later arrival 
j > i in the following steps. As a result, every interval will yield at most one false 
computation. Since there are exactly N intervals this completes the proof. |j 

Proof of Lemma 15.21 Since x > ai it holds that i < 7i ri and thus using (15. 2ft if 
ki < 7 r r then 6i = 0, and if ki > n r then 

/3 — a (hi — i ) 

Since N < /d/a + 1, the denominator is always positive. We next show that the 
numerator is non-negative by induction on h 7Tr < i < ix r . Recall that hi = min{h : 
kh > i}, and so ki > n r is equivalent to i > h nr . Thus for j < h 7Tr : 9j = 0 and 
the denominator in the case i = h 7Tr equals a(l — 0) > 0. The induction step is then 
immediate because the sum ^ = / t . 9j is non-negative for all h nr < i < 7r r . | 

Proof of Lemma 15.31 Without loss of generality assume that a + ^ < a* < a — Vi E 
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A f. If this were not the case we could always extend the search range by in both 
directions. Hence, Qi = 0 at x = a and at x = a for any i e n. Furthermore, from 
Lemma [5.21 we have that < 0 for x > dj. Clearly, there is some x such that > 0 
for the first time. So far we have established that 6 t starts at zero, is positive at some 
point and negative at some back to zero, for every i 6 n. We are left with finding the 
number of possible sign changes prior to x — a % . For any x < a,i it follows that i > n r , 
and from (15. 2 p we have that: 


Oi 


(3 — a(ki — i ) 


Note that 6 t can only be negative when there is at least one j < i such that 9j < 0. 
We use this to complete the proof by induction on the initial order tt. We start the 
induction at i = 2 because n r = 1 in the initial permutation and d 7Tr > 0 for all values 
of x. For i — 2 and x < <i 2 '. O 2 = a ^a{k 2 ~X) — 0- Together with Lemma [5721 we have 
established that 62 changes sign exactly once. Now let us assume that the claim is 
correct for all j < i — 1. From (15. 2\i we see that for x < a*, 6 j can only change sign 
when one of the previous j G {h il ..., i — 1} changes sign. If 6 W 1 changed sign exactly 
i — 2 times then 9i can potentially change at all these times and additionally when 
x = a*, and therefore there are indeed at most i — 1 changes of sign. | 
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